1 Exercicis

Els exercicis es realitzaran sobre la base del joc de dades Hawks present en el paquet R Stat2Data.

Els estudiants i el professorat del Cornell College a Mount Vernon, Iowa, van recollir dades durant molts anys al mirador de falcons de l’estany MacBride, prop d’Iowa City, a l’estat d’Iowa. El joc de dades que analitzem aquí és un subconjunt del conjunt de dades original, utilitzant només aquelles espècies per a les que hi havia més de 10 observacions. Les dades es van recollir en mostres aleatòries de tres espècies diferents de falcons: Cua-roja, Esparver i Falcó de Cooper.

Hem seleccionat aquest joc de dades per la seva semblança amb el joc de dades penguins i pel seu potencial alhora d’aplicar-li algoritmes de mineria de dades no supervisats. Les variables numèriques en què us basareu són: Wing, Weight, culmen, Hallux

if (!require('Stat2Data')) install.packages('Stat2Data'); library('Stat2Data')
data("Hawks")
summary(Hawks)
##      Month             Day             Year       CaptureTime   ReleaseTime 
##  Min.   : 8.000   Min.   : 1.00   Min.   :1992   11:35  : 14          :842  
##  1st Qu.: 9.000   1st Qu.: 9.00   1st Qu.:1995   13:30  : 14   11:00  :  2  
##  Median :10.000   Median :16.00   Median :1999   11:45  : 13   11:35  :  2  
##  Mean   : 9.843   Mean   :15.74   Mean   :1998   12:10  : 13   12:05  :  2  
##  3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.:2001   14:00  : 13   12:50  :  2  
##  Max.   :11.000   Max.   :31.00   Max.   :2003   13:05  : 12   13:32  :  2  
##                                                  (Other):829   (Other): 56  
##       BandNumber  Species  Age     Sex          Wing           Weight      
##            :  2   CH: 70   A:224    :576   Min.   : 37.2   Min.   :  56.0  
##  1142-09240:  1   RT:577   I:684   F:174   1st Qu.:202.0   1st Qu.: 185.0  
##  1142-09241:  1   SS:261           M:158   Median :370.0   Median : 970.0  
##  1142-09242:  1                            Mean   :315.6   Mean   : 772.1  
##  1142-18229:  1                            3rd Qu.:390.0   3rd Qu.:1120.0  
##  1142-19209:  1                            Max.   :480.0   Max.   :2030.0  
##  (Other)   :901                            NA's   :1       NA's   :10      
##      Culmen         Hallux            Tail        StandardTail  
##  Min.   : 8.6   Min.   :  9.50   Min.   :119.0   Min.   :115.0  
##  1st Qu.:12.8   1st Qu.: 15.10   1st Qu.:160.0   1st Qu.:162.0  
##  Median :25.5   Median : 29.40   Median :214.0   Median :215.0  
##  Mean   :21.8   Mean   : 26.41   Mean   :198.8   Mean   :199.2  
##  3rd Qu.:27.3   3rd Qu.: 31.40   3rd Qu.:225.0   3rd Qu.:226.0  
##  Max.   :39.2   Max.   :341.40   Max.   :288.0   Max.   :335.0  
##  NA's   :7      NA's   :6                        NA's   :337    
##      Tarsus        WingPitFat        KeelFat           Crop       
##  Min.   :24.70   Min.   :0.0000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:55.60   1st Qu.:0.0000   1st Qu.:2.000   1st Qu.:0.0000  
##  Median :79.30   Median :1.0000   Median :2.000   Median :0.0000  
##  Mean   :71.95   Mean   :0.7922   Mean   :2.184   Mean   :0.2345  
##  3rd Qu.:87.00   3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:0.2500  
##  Max.   :94.00   Max.   :3.0000   Max.   :4.000   Max.   :5.0000  
##  NA's   :833     NA's   :831      NA's   :341     NA's   :343

1.1 Exercici 1

Presenta el joc de dades, nom i significat de cada columna, així com les distribucions dels seus valors.

Adicionalment realitza un estudi similar al dels exemples 1.1 i 1.2

1.1.1 Exploració de les dades


El joc de dades ‘Hawks’ conté 908 registres de 19 variables. Les variables corresponen a les següents observacions:

  • Month: mes de l’observació.

  • Day: dia.

  • Year: any.

  • CaptureTime: hora de la captura (HH:MM)

  • ReleaseTime: hora de l’alliberament (HH:MM)

  • BandNumber: codi d’identitat a l’anell identificador.

  • Species: CH=Falcó de Cooper, RT=Cua-roja, SS=Esparver

  • Age: A=Adult or I=Immadur

  • Sex: F=Femella or M=Mascle

  • Wing: llargada en mm. de la ploma principal de l’ala, des de la punta fins l’articulació.

  • Weight: pes de l’ocell en gr.

  • Culmen: llargada en mm. de la cresta externa del bec.

  • Hallux: llargada en mm. del dit posterior.

  • Tail: mesurament propi de la llargada de la cua en mm.

  • StandardTail: mesurament estàndard de la llargada de la cua en mm.

  • Tarsus: llargada de l’os principal del peu en mm.

  • WingPitFat: quantitat de greix a l’aixella

  • KeelFat: quantitat de greix al pit

  • Crop: quantitat de material al pap. 1=ple, 0=buit.


Primer de tot, analitzarem visualment les distribucions de les 4 variables d’interès per detectar la presència de valor aïllats:

boxplot(Hawks$Wing, main = 'Wing')

boxplot(Hawks$Weight, main = 'Weight')

boxplot(Hawks$Culmen, main = 'Culmen')


Veiem que a les variables ‘Wing’, ‘Weight’ i ‘Culmen’ no hi ha valors aïllats. A la variable ‘Hallux’ però, sí que hi veiem un parell de valors molt allunyats:

boxplot(Hawks$Hallux, main = 'Hallux')


Examinem aquests dos valors:

na.omit(Hawks[Hawks$Hallux > 300, 10:13])
##    Wing Weight Culmen Hallux
## 91  410   1210   27.5  308.0
## 96  418   1180   30.1  341.4

Veiem que en aquests dos casos la resta de variables també presenta valors grans, situats en el quartil superior. Per tant, assumim que els valors no són erronis, i que es tracta simplement d’exemplars d’animals de grans dimensions.


Seguidament, explorarem la presència de valors nuls:

nrow(Hawks[is.na(Hawks$Weight) | is.na(Hawks$Wing) | is.na(Hawks$Culmen) | is.na(Hawks$Hallux),])
## [1] 17
nrow(Hawks[is.na(Hawks$Weight) | is.na(Hawks$Wing) | is.na(Hawks$Culmen) | is.na(Hawks$Hallux),]) / nrow(Hawks)
## [1] 0.01872247


Malgrat que la quantitat de files amb valors nuls en les nostres variables d’interès només representen un 1.8% del total de registres, creiem que seria interessant no deixar perdre aquestes registres, i plantegem imputar els valors amb el mètode de veïns més propers:

library(VIM)
## Loading required package: colorspace
## Loading required package: grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep


Fem la imputació dels valors, utilitzant les mateixes variables d’interès com a referència:

hawks_knn_imp <- kNN(Hawks, variable = c("Weight", "Wing", "Culmen", "Hallux"), dist_var = c("Weight", "Wing", "Culmen", "Hallux"), imp_var = FALSE)


Creem una variable amb els registres que contenien valors nuls per a qualsevol de les nostres variables:

hawks_na <- which(is.na(Hawks$Weight) | is.na(Hawks$Wing) | is.na(Hawks$Culmen) | is.na(Hawks$Hallux))


Mostrem una taula amb els valors abans i després de la imputació:

library(knitr)
hawks_table_na <- Hawks[hawks_na, c("Weight", "Wing", "Culmen", "Hallux")]
hawks_table_imp <- hawks_knn_imp[hawks_na, c("Weight", "Wing", "Culmen", "Hallux")]
colnames(hawks_table_imp)[colnames(hawks_table_imp) %in% c("Weight", "Wing", "Culmen", "Hallux")] <- c("Weight_imp", "Wing_imp", "Culmen_imp", "Hallux_imp")
hawks_imputed <- cbind(hawks_table_na, hawks_table_imp)
hawks_imputed <- hawks_imputed[, c("Weight", "Weight_imp", "Wing", "Wing_imp","Culmen", "Culmen_imp", "Hallux", "Hallux_imp")]
kable(hawks_imputed)
Weight Weight_imp Wing Wing_imp Culmen Culmen_imp Hallux Hallux_imp
2 930 930 376 376 NA 26.0 NA 30.8
14 NA 1170 393 393 28.2 28.2 30.6 30.6
26 100 100 173 173 NA 11.4 NA 12.2
64 NA 971 326 326 25.2 25.2 27.7 27.7
69 NA 170 194 194 11.4 11.4 14.2 14.2
194 1040 1040 403 403 NA 26.1 29.9 29.9
218 NA 945 202 202 NA 24.5 NA 31.2
263 480 480 NA 261 17.7 17.7 32.1 32.1
315 NA 945 162 162 NA 24.5 NA 31.2
332 1244 1244 368 368 NA 27.5 NA 30.9
407 NA 895 361 361 24.4 24.4 27.9 27.9
414 NA 1055 271 271 27.4 27.4 33.0 33.0
420 NA 185 205 205 13.7 13.7 15.0 15.0
433 1080 1080 381 381 NA 26.3 32.3 32.3
523 NA 170 190 190 11.9 11.9 14.6 14.6
525 NA 1240 406 406 27.0 27.0 31.1 31.1
886 945 945 363 363 24.5 24.5 NA 27.5


Finalment, mostrem el sumari del dataset amb el qual treballarem, amb 908 observacions i 4 variables:

hawks4 <- hawks_knn_imp[, 10:13]
summary(hawks4)
##       Wing           Weight           Culmen          Hallux      
##  Min.   : 37.2   Min.   :  56.0   Min.   : 8.60   Min.   :  9.50  
##  1st Qu.:202.0   1st Qu.: 185.0   1st Qu.:12.80   1st Qu.: 15.10  
##  Median :370.0   Median : 970.0   Median :25.50   Median : 29.40  
##  Mean   :315.6   Mean   : 772.1   Mean   :21.82   Mean   : 26.42  
##  3rd Qu.:390.0   3rd Qu.:1120.0   3rd Qu.:27.30   3rd Qu.: 31.40  
##  Max.   :480.0   Max.   :2030.0   Max.   :39.20   Max.   :341.40


1.1.2 Agregació amb k-means


Instal·lem la llibreria:

if (!require('cluster')) install.packages('cluster')
## Loading required package: cluster
library(cluster)


Comencem a cercar el nombre òptim de clústers. Provem primer amb la funció kmeans, iterant per agrupacions d’entre 2 i 10 clústers:

d <- daisy(hawks4) 
results <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
  fit           <- kmeans(hawks4, i)
  y_cluster     <- fit$cluster
  sk            <- silhouette(y_cluster, d)
  results[i] <- mean(sk[,3])
}
library(NbClust)
library(factoextra)
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_nbclust(hawks4, pam, method = "silhouette")+ theme_classic()

En el gràfic podem veure que l’índex de silhouette dóna com a òptima l’agregació de 2 clústers. Mostrem a continuació el detall d’aquests índexs:


fit2       <- kmeans(hawks4, 2)
y_cluster2 <- fit2$cluster

fit3       <- kmeans(hawks4, 3)
y_cluster3 <- fit3$cluster

fit4       <- kmeans(hawks4, 4)
y_cluster4 <- fit4$cluster

fit5       <- kmeans(hawks4, 5)
y_cluster5 <- fit5$cluster
sk2 <- silhouette(y_cluster2, d)
sk3 <- silhouette(y_cluster3, d)
sk4 <- silhouette(y_cluster4, d)
sk5 <- silhouette(y_cluster5, d)

mean(sk2[,3])
## [1] 0.8057608
mean(sk3[,3])
## [1] 0.6623463
mean(sk4[,3])
## [1] 0.6371349
mean(sk5[,3])
## [1] 0.5945413


Seguidament, com a alternativa, generem un gràfic amb el càlcul de la suma de quadrats de distàncies respecte a centroides, també per a entre 2 i 10 agrupacions:


results <- rep(0, 10)
for (i in c(2,3,4,5,6,7,8,9,10))
{
  fit           <- kmeans(hawks4, i)
  results[i] <- fit$tot.withinss
}
plot(2:10,results[2:10],type="o",col="magenta",pch=0,xlab="Number of clusters",ylab="tot.tot.withinss")

Aquest gràfic resulta més difícil d’interpretar, ja que no hi trobem un canvi d’inclinació molt pronunciat abans de k=6, i a partir de k=6 hi ha pujades i baixades.

Si seguim explorant alternatives, amb els índexs ‘Calinski-Harabasz’ (‘ch’) i silueta mitjana (‘asw’) obtenim una quantitat de clústers òptima de 10 i 2 clústers respectivament:


if (!require('fpc')) install.packages('fpc'); library('fpc')
## Loading required package: fpc
fit_ch  <- kmeansruns(hawks4, krange = 1:10, criterion = "ch") 
fit_asw <- kmeansruns(hawks4, krange = 1:10, criterion = "asw") 
fit_ch$bestk
## [1] 10
fit_asw$bestk
## [1] 2


plot(1:10,fit_ch$crit,type="o",col="green",pch=0,xlab="Number of clusters",ylab="Calinski-Harabasz")


plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Number of clusters",ylab="mean silhouette")


Havent fet totes aquestes simulacions, veiem que l’algoritme no acaba d’identificar com a òptim el nombre de clústers (3) que nosaltres sabem que equivaldria a la divisió en espècies de la mostra. En la majoria dels casos, els índexs mostren com a òptima una agrupació en 2 clústers.


1.1.3 Resultats comparats


Aprofitant doncs que sabem que la mostra es divideix en 3 classes (espècies) diferents, ara visualitzarem la comparativa entre agrupacions fetes amb k-means (k=3) i les dades agrupades segons espècie.


Primer de tot visualitzem la relació entre les variables ‘Wing’ i ‘Weight’:

hawks3clusters <- kmeans(hawks4, 3)

plot(hawks4[c(1,2)], col=hawks3clusters$cluster, main="Classificació k-means")


hawks5 <- hawks_knn_imp[, c(7,10:13)]

plot(hawks4[c(1,2)], col=as.factor(hawks5$Species), main="Classificació real")

En aquesta comparació veiem l’agrupació que fa l’algoritme no es correspon de forma molt acurada a la divisió per espècies. L’espècie més nombrosa (marcada en vermell al segon gràfic) queda dividida en dos grups per l’algoritme k-means, mentre que les altres dues espècies s’agrupen en un sol clúster.

Farem la mateixa comparació per a diversos parells de variables, i trobarem en tots els casos el mateix tipus de discrepància entre l’algoritme i la classificació real: l’algoritme tendeix a dividir l’espècie ‘RT’ en dos grups, i ajunta les altres dues espècies en un sol grup.


plot(hawks4[c(3,4)], col=hawks3clusters$cluster, main="Classificació k-means")


plot(hawks4[c(3,4)], col=as.factor(hawks5$Species), main="Classificació real")


plot(hawks4[c(1,3)], col=hawks3clusters$cluster, main="Classificació k-means")

plot(hawks4[c(1,3)], col=as.factor(hawks5$Species), main="Classificació real")


plot(hawks4[c(1,4)], col=hawks3clusters$cluster, main="Classificació k-means")

plot(hawks4[c(1,4)], col=as.factor(hawks5$Species), main="Classificació real")


plot(hawks4[c(2,3)], col=hawks3clusters$cluster, main="Classificació k-means")

plot(hawks4[c(2,3)], col=as.factor(hawks5$Species), main="Classificació real")


plot(hawks4[c(2,4)], col=hawks3clusters$cluster, main="Classificació k-means")

plot(hawks4[c(2,4)], col=as.factor(hawks5$Species), main="Classificació real")




1.2 Exercici 2


1.2.1 Generació del model DBSCAN


Carreguem les llibreries necessàries:

if (!require('dbscan')) install.packages('dbscan'); library('dbscan')
## Loading required package: dbscan
## 
## Attaching package: 'dbscan'
## The following object is masked from 'package:fpc':
## 
##     dbscan
## The following object is masked from 'package:VIM':
## 
##     kNN


Per tal d’utilitzar l’algorisme DBSCAN, primer de tot intentarem trobar la mida òptima del radi (valor de eps) per la nostra mostra.


Primer de tot aplicarem l’algorisme amb un valor arbitrari (5 en aquest cas):

tryeps <- dbscan(hawks4, eps = 5)
tryeps
## DBSCAN clustering for 908 objects.
## Parameters: eps = 5, minPts = 5
## The clustering contains 12 cluster(s) and 663 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12 
## 663  95  66   6   9   6   5  15   5  16   7   9   6 
## 
## Available fields: cluster, eps, minPts

Veiem que l’algorisme genera 11 clústers i 663 valors allunyats, i per tant no en resulta una classificació gens idònia.


Per tal de trobar un valor d’eps que ens doni resultats útils, generarem un gràfic de distància k amb les nostres dades, i identificarem el punt d’inflexió del gràfic per trobar un valor òptim de eps:

kNNdistplot(hawks4, k = 5)
abline(h=50, col = "red", lty=2)

En el gràfic hem superposat una línia vermella en el valor 50, que ens sembla aproximadament el valor on hi ha una inflexió. Utilitzarem doncs el valor eps=50 per implementar l’algorisme. El valor minPts el deixarem en 5 (4 variables + 1):

opt <- dbscan(hawks4, eps = 50, MinPts = 5)
## Warning in dbscan(hawks4, eps = 50, MinPts = 5): converting argument MinPts
## (fpc) to minPts (dbscan)!
opt
## DBSCAN clustering for 908 objects.
## Parameters: eps = 50, minPts = 5
## The clustering contains 3 cluster(s) and 28 noise points.
## 
##   0   1   2   3 
##  28 561  64 255 
## 
## Available fields: cluster, eps, minPts

Veiem com l’algorime identifica 3 grups, i determina que hi ha 28 observacions no classificables (soroll)


Ara apliquem l’algorisme optics per tal de generar una variable que utilitzarem més avall. Veiem que el resultat d’extreure’n un dbscan és el mateix:

opti <- optics(hawks4)
dbs <- extractDBSCAN(opti, eps_cl = 50)
dbs
## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 5, eps = 446.7443564277, eps_cl = 50, xi = NA
## The clustering contains 3 cluster(s) and 28 noise points.
## 
##   0   1   2   3 
##  28 561  64 255 
## 
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
##                   eps_cl, xi, cluster


Ara visualitzem el resultat. Els tres colors a les valls representen els diferents grups, i els pics de color negre són els valors outliers:

plot(dbs)


Ara, fem un gràfic que mostra quins són els clústers que proposa l’algorisme dbscan quan visualitzem les observacions per parelles de variables:

pairs(hawks5, col= opt$cluster)


Per tal de comprovar si aquestes agrupacions s’ajusten a la classificació real, generem el mateix tipus de gràfic utilitzant els colors per indicar les espècies. Comparant els dos gràfics, veiem que les agrupacions són força similars:

pairs(hawks5, col= hawks5$Species)



Per veure-ho més en detall, mostrem un gràfic amb només dues variables, fent la comparativa entre la classificació DBSCAN i la classificació real:

plot(hawks4[c(1,2)], col=opt$cluster, main="Classificació DBSCAN")

plot(hawks4[c(1,2)], col= as.factor(hawks5$Species), main="Classificació real")


Veiem que, a grans trets, la classificació que fa DBSCAN és molt propera a la classificació real. La diferència més notable en el gràfic és que DBSCAN no inclou els valors considerats com a outliers, que es troben fora del radi definit i que per tant resulten inclassificables (també hi ha alguns valors que queden allunyats del centre del seu grup que es classifiquen erròniament en un grup diferent). Generem un altre gràfic on es mostren aquests valors que queden fora del radi (punts en negre):

hullplot(hawks4[c(1,2)], dbs)



Presentem una última comparativa entre el model i la classificació real, aquest cop per les variables ‘weight’ i ‘culmen’, i amb la llibreria ggplot, on podem veure les etiquetes de cada espècie en el gràfic de la classificació real, i també els punts aïllats (etiqueta ‘0’) al gràfic del model DBSCAN:

library(ggplot2)
qplot(Weight, Culmen, data = hawks5, colour = Species, main = 'Classificació Real')


qplot(Weight, Culmen, data = hawks4, colour = as.factor(opt$cluster), main = "Classificació DBSCAN")


1.2.2 Avaluació del model



Finalment, realitzarem una avaluació del model DBSCAN que hem utilitzat. Per fer-ho, utilitzarem l’índex de silhouette aplicat sobre la comparació entre el nostre model i una matriu de distàncies de les dades originals de la qual haurem extret els valors que el model havia identificat com a outliers:

noise <- opt$cluster==0
clusters <- opt$cluster[!noise]
d <- dist(hawks4[!noise, 1:4])
silh <- silhouette(clusters, d)
plot(silh, border=NA, col=sort(clusters), main="")



Veiem en el resultat que la mitjana de l’índex de silhouette és de 0.74, un índex que denota una bona qualitat d’agrupament (valors de silhouette entre -1.0 i 1.0).





1.3 Exercici 3

Realitza una comparativa dels mètodes k-means i DBSCAN

1.3.1 Resposta 3


D’una banda, els resultats de la classificació amb l’algorisme k-means obtenien una avaluació interna òptima quan la mostra es classificava en 2 clústers. Tanmateix, sabem que la mostra es divideix en tres grups (espècies), i per tant hem comprovat que k-means no és el mètode òptim per revelar l’estructura de les nostres dades. Concretament, quan forçàvem l’algorisme k-means a formar 3 clústers, aquest tendia a dividir l’espècie més nombrosa en dos clústers diferents, i agrupava les altres dues espècies en una de sola. K-means és particularment sensible als valors aïllats, i possiblement això creava aquesta distorsió en relació a la classificació real.

Aquesta sensibilitat respecte als valors aïllats s’evita amb l’algorisme DBSCAN, però la contrapartida és que els valors aïllats es classifiquen com a soroll i queden descartats del model. Això és degut al fet que l’algorisme basa la seva classificació en la densitat de les agrupacions, i no pot incorporar de forma eficient els valors que escapen als radis de proximitat que són necessaris per definir els clústers. En el nostre cas hem aconseguit minimitzar el nombre de valors exclosos a 28 valors sobre un total de 908; és una pèrdua petita, però no del tot negligible.

Pel que fa a la comparativa entre els dos mètodes, amb DBSCAN cal fer el treball extra de recercar els valors òptims de radi de veïnatge, mentre que amb k-means cal definir d’entrada el nombre de clústers en el qual volem classificar la mostra. Amb les eines que hem utilitzat per definir aquests factors per ambdós algorismes (diferents mesures d’avaluació interna per a k-means, que donaven com a òptima una divisió en 2 clùsters, i l’anàlisi visual d’un gràfic de distàncies k per a DBSCAN, que ens donava un valor òptim d’eps entorn a 50), i comparant els resultats dels algorismes amb la classificació real, podem afirmar que, almenys en aquest cas, l’algorisme DBSCAN ha estat més efectiu a l’hora d’apropar-se a la classificació real de les dades.